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Abstract 

One  major  issue  in  the  accurate  solution  of  convection-dominated  problems 
by  means  of  high-order  methods  is  the  ability  of  the  solver  to  maintain 
monotonicity.  This  problem  is  critical  for  spectral  elements,  where  Gibbs 
oscillations  may  pollute  the  solution.  However,  typical  filter-based  stabi¬ 
lization  techniques  used  with  spectral  elements  are  not  monotone.  In  this 
paper,  residual-based  stabilization  methods  originally  derived  for  finite  ele¬ 
ments  are  constructed  and  applied  to  high-order  spectral  elements.  In  par¬ 
ticular,  we  show  that  the  use  of  the  Variational  Multiscale  (VMS)  method 
greatly  improves  the  solution  of  the  transport-diffusion  equation  by  reduc¬ 
ing  over-  and  under-shoots,  and  can  be  therefore  considered  an  alternative 
to  the  limitations  of  filter-based  schemes.  We  also  combine  these  methods 
with  discontinuity  capturing  schemes  to  suppress  oscillations  that  may  oc¬ 
cur  in  proximity  of  boundary  or  internal  layers.  Additional  improvement  in 
the  solution  is  also  obtained  when  p-adaptivity  is  used  in  combination  with 
VMS  in  the  regions  where  discontinuities  occur.  The  algorithms  are  assessed 
with  the  solution  of  classical  steady  and  transient  one-  and  two-dimensional 
problems  using  spectral  elements  up  to  order  16. 
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1.  Introduction 

A  large  number  of  physical  applications  relies  on  the  accurate  solution 
of  the  transport-diffusion  equation 

£(q)  ■=  +  u  '  Vg  -  vAq  =  /,  (1) 

where  v  >  0  is  a  diffusion  coefficient,  u  is  a  known  velocity  field,  and  the  un¬ 
known  q  is  the  tracer  to  be  transported.  The  solution  of  (1)  should  respect 
two  significant  properties:  (i)  positivity  must  be  preserved,  and  (ii)  smearing 
at  internal  and  boundary  layers  should  not  be  excessive.  These  properties 
are  extremely  important  in  the  context  of  transport  in  the  atmosphere.  Both 
limited-area  and  global  atmospheric  models  for  weather  prediction  need 
monotonic  advection  of  tracers  and  moisture  variables,  otherwise  the  wrong 
amount  of  precipitation  would  be  forecasted.  Simple  nricrophysics  schemes, 
such  as  Kessler  [1],  require  three  variables  (water  vapor,  cloud  water,  and 
rain) ,  whereas  more  sophisticated  parameterizations  include  additional  vari¬ 
ables  such  as  ice  and  snow  [2] .  Similarly,  climate  models  require  transport  of 
hundreds  of  tracers,  each  representing  a  different  chemical  species.  Regard¬ 
less  of  the  physical  scales  of  the  model,  tracers  must  remain  positive  since 
the  physical  parameterizations  that  govern  sub-grid  scale  processes  such 
as  auto-conversion  and  sedimentation,  implicitly  assume  such  a  condition. 
These  issues  have  been  addressed  for  both  transient  and  stationary  prob¬ 
lems  (See,  e.g.,  [3])  and,  in  the  context  of  finite  element  methods,  so-called 
stabilized  methods  have  been  an  active  topic  of  research  since  their  intro¬ 
duction  in  the  early  1980s  with  the  Streamline-Upwind  method  of  Hughes 
and  Brooks  [4].  In  this  paper  we  address  the  problem  of  solving  (1)  by 
high-order  spectral  element  methods  without  losing  the  ability  to  approach 
a  monotone  solution  to  the  problem.  Higher-order  accuracy,  in  fact,  comes 
at  the  price  of  aliasing  phenomena  in  the  solution  ( [5] ) ,  but  the  anti-aliasing 
filters  typically  used  to  give  a  stable  spectral  element  solution  do  not  respect 
the  two  conditions  described  above.  Therefore,  to  achieve  monotonic  results 
with  high-order  spectral  elements,  we  consider  stabilization  methods  origi¬ 
nally  devised  for  finite  elements,  and  focus  on  methods  that  can  be  derived 
directly  from  subgrid  scale  considerations  as  originally  defined  in  [6]  and  [7] 
in  the  context  of  variational  multiscale  methods.  These  techniques  assure 
stability  by  designing  a  diffusion-type  term  that  is  added  to  the  Galerkin 
formulation  of  the  original  problem. 

The  first  stabilized  schemes  based  on  the  addition  of  a  diffusive  stabi¬ 
lization  term  to  the  Galerkin  equation  are  the  Artificial  Viscosity  methods 
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(AV)  [8]  and  the  Streamline-  Upwind  method  (SU)  [4].  AV,  as  Hyper-viscosity 
(HV),  is  often  used  in  atmospheric  and  ocean  modeling  due  to  the  property 
of  preserving  the  correct  energy  cascade  in  simulations  that  involve  turbu¬ 
lence.  The  SU  scheme  uses  the  information  in  the  direction  of  the  flow  to 
add  viscosity  only  in  the  streamline  direction.  Both  methods  use  a  constant 
diffusion  coefficient  that  does  not  typically  change  from  element  to  element. 
A  major  improvement  came  by  introducing  the  residual  of  the  governing 
equation  in  the  definition  of  the  stabilization  term.  When  the  computed 
solution  approaches  the  exact  solution,  the  stabilization  term  should  vanish. 
These  schemes,  which  are  consistent  in  that  the  stabilization  terms  goes  to 
zero  as  the  numerical  solution  approaches,  are  considered  in  this  paper.  The 
most  commonly  used  are  the  Upwind- Streamline/ Petrov- Galerkin  (SUPG) 
and  the  Galerkin/ Least- Squares  (GLS),  devised  in  1982  [9]  and  1989  [10], 
respectively,  as  a  consistent  counterpart  to  SU.  GLS  was  designed  as  a  gener¬ 
alization  of  SUPG,  but  in  the  limit  of  pure  advection,  or  for  piece- wise  linear 
elements,  the  GLS  and  SUPG  methods  are  equivalent.  Stability  analysis  for 
these  two  methods  is  detailed  in  [11,  12,  10].  The  Gradient  Galerkin/ Least- 
Squares  [13]  for  advection-diffusion  with  a  reaction  term,  or  the  Unusual 
Stabilized  Finite  Element  Method  (USFEM)  [14,  15]  are  a  few  examples.  In 
the  framework  of  high  order  methods,  Petrov-Galerkin  stabilization  was  ap¬ 
plied  by  Pasquarelli  and  Quarteroni  [16]  to  stabilize  the  convection-diffusion 
equation  with  the  spectral  method.  Canuto  used  bubble  functions  to  address 
the  same  issue  [17]  (See  also  [18,  19]). 

The  analyses  of  Hughes  [6] ,  Hughes  and  Stewart  [20] ,  and  Hughes  et  al. 
[7]  form  the  unifying  theory  of  all  stabilized  finite  element  methods.  Sta¬ 
bilized  methods  are  subgrid  scale  models  where  the  unresolved  scales  are 
intimately  related  to  the  instabilities  at  the  level  of  the  resolved  scales,  and 
thus  should  be  used  in  the  construction  of  the  stabilization  term.  These 
schemes  are  known  as  Variational  Multiscale  (VMS)  methods.  Details  are 
given  in  subsection  2.2.2.  VMS  methods  are  all  residual-based  methods  that 
improve  the  stability  properties  of  the  solution,  and  preserve  the  accuracy 
of  the  underlying  numerical  scheme  [21],  However,  Godunov’s  theorem  [22] 
implies  that  the  latter  property  may  be  violated  in  the  proximity  of  discon¬ 
tinuities  or  strong  gradients.  As  already  observed  by  Hughes,  Cottrell,  and 
Bazilevs  in  [21],  where  NURBS  were  used  as  high-order  basis  functions,  un¬ 
expected  convergence  to  monotone  results  were  obtained  independently  of 
the  order  of  the  polynomials  used.  Neither  SUPG,  GLS,  nor  VMS,  however, 
preclude  the  formation  of  over-  and  under-shoots  in  the  proximity  of  sharp 
gradients  of  the  solution.  For  this  reason,  shock  capturing  or  discontinuity 
dissipation  techniques,  also  referred  to  as  Spurious  oscillations  at  layers  di- 
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minishing  methods  (SOLD)  are  used  in  combination  with  SUPG  and  VMS 
to  introduce  an  additional  term  to  the  stabilized  form  of  the  equation.  This 
issue  was  treated  for  the  first  time  in  [23],  where  details  on  how  to  build 
the  stabilization  parameter  are  also  given,  and  in  [24]  for  non-linear  prob¬ 
lems.  A  detailed  review  of  most  existing  SOLD  schemes  can  be  found  in 
a  two-part  paper  by  John  and  Knobloch  [25,  26],  where  a  modification  of 
the  discontinuity -capturing  of  Codina  [27]  is  presented  and  is  shown  to  be  a 
promising  option  for  FE  solutions  characterized  by  boundary  layers. 

All  these  methods  strongly  depend  on  a  stabilization  parameter  that  will 
be  identified  by  r  throughout  the  paper.  It  will  be  also  referred  to  as  in¬ 
trinsic  time.  A  classical  result  was  obtained  by  Franca,  Frey  and  Hughes  in 
[28]  by  error  analysis.  Their  result  was  reproduced  by  other  authors  using 
different  approaches.  Additional  expressions  for  r  were  found  by  Codina  in 
[29,  30],  by  Codina,  Onate  and  Cervera  in  [31],  by  Harari  and  Hughes  in 
[13],  and  by  Shakib,  Hughes  and  Johan  in  [32],  who  based  the  derivation 
on  the  (discrete)  maximum  principle.  Another  expression  is  due  to  Franca 
and  Valentin  [15]  who  based  their  derivation  on  convergence  and  stability 
analysis.  Starting  with  the  formalization  of  VMS  methods  by  Hughes  [6] ,  r 
has  often  been  derived  using  Green’s  functions,  a  thorough  analysis  of  which 
is  done  by  Hughes  and  Sangalli  in  [33] .  Recently,  Houzeaux,  Eguzkitza  and 
Vazquez  [34]  proposed  a  new  way  to  derive  the  approximate  subgrid  scale 
solution,  with  results  that  are  comparable  to  those  of  Hauke  and  Garcia- 
Olivares  in  [35].  In  [36],  Codina  builds  t  using  the  Fourier  analysis  of  the 
problem;  however,  the  correct  choice  of  r  remains  an  open  problem.  For  this 
reason,  we  propose  r  for  higher-order  spectral  elements  and  use  it  to  con¬ 
struct  an  appropriate  stabilization  method.  To  further  improve  the  solution, 
we  combine  VMS  stabilization  with  a  p-adaptivity  algorithm,  p-adaptivity 
(see,  e.g.,  [37,  38])  adjusts  the  order  of  the  polynomial  interpolation  instead 
of  modifying  the  computational  grid.  In  this  paper  we  use  it  in  combina¬ 
tion  with  VMS  in  those  regions  where  discontinuities  cause  the  high-order 
solution  to  be  affected  by  Gibbs  oscillations. 

1.1.  Main  contribution  of  this  paper 

The  main  problem  that  we  want  to  solve  with  the  work  presented  here 
is  that  of  stabilizing  the  spectral  element  solution  of  advection-dominated 
problems  by  sub-grid  scale  stabilization  techniques  (namely,  VMS),  and  im¬ 
prove  the  solution  by  reducing  the  under-  and  overshoots  that  would  occur 
if  classical  filters  were  to  be  used.  The  definition  and  implementation  of  r 
in  VMS  stabilization  may  greatly  affect  the  result.  We  hence  adapted  the 
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method  described  in  [34]  to  compute  r,  and  applied  it  to  spectral  elements 
that  use  Legendre-Gauss-Lobatto  (LGL)  nodes,  and  show  that  this  tech¬ 
nique  can  be  used  in  problems  where  spectral  element  filters  fail.  We  finally 
apply  polynomial  adaptivity  (i.e.,  we  lower  the  order  of  interpolation  by 
keeping  the  computational  grid  untouched) ,  only  where  the  solution  is  char¬ 
acterized  by  a  propagating  discontinuity.  The  combination  of  VMS,  shock 
capturing  and  p-adaptive  methods  will  be  shown  to  be  an  encouraging  di¬ 
rection  to  take  for  constructing  high-order  positive-definite  spectral  element 
methods.  To  assess  the  algorithm,  steady  and  transient  advection-diffusion 
problems  are  solved  on  one-  and  two-dimensional  domains  using  spectral 
elements  up  to  order  1 6.  We  compare  the  performance  of  the  method  using 
the  VMS  against  those  obtained  with  previous  classical  SE  schemes. 

1.2.  Outline  of  this  paper 

The  remainder  of  the  paper  is  organized  as  follows.  After  this  introduc¬ 
tion,  the  numerical  method  and  the  corresponding  stabilization  are  derived 
in  Sections  2  and  3.  Tests  to  verify  the  algorithm  and  a  discussion  of  the 
results  are  presented  in  Sections  4  and  5,  respectively. 

2.  Numerical  method 

Given  the  space  L 2  of  real-valued  functions  that  are  square  integrable 
in  a  bounded  domain  fl  C  l2  with  boundary  T,  the  Sobolev  space  H 1  of 
weakly-differentiable  functions  will  be  used.  Specifically,  W  C  H 1  represents 
the  space  of  trial  and  basis  functions  of  the  Galerkin  formulation  to  follow. 
In  L2  the  inner  product  is  given  by  (•,  •),  and  the  ^-norm  associated  with 
the  space  is  denoted  by  ||  •  || 2 -  For  simplicity,  we  add  the  property  that  the 
solution  q  vanishes  on  the  boundary  <917;  under  this  assumption,  W  C  Hq. 
Given  a  finite  element  partition  19/,  =  U^i  Tv,  of  the  computational  domain 
19  into  nei  high-order  conforming  quadrilaterals  of  characteristic  length  h, 
Wh  is  the  finite  dimensional  set  derived  from  W.  The  discrete  weak  form 
reduces  to  the  problem  of  finding  the  function  qh  €  [Wh\  0,  t)  such  that 

(^,^)  +  a(^,/)  =  (^,/)  y^hewh  (2) 

where,  after  integrating  by  parts  and  assuming  homogeneous  Dirichlet  bound¬ 
ary  conditions,  we  define: 


8qh 
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a{i/>h,qh) 


ip hu  ■  X7qh  dVth  + 


vViph  •  \7qhdnh, 


(Tph,f)=  [  *Phfdnh  =  o. 

JQh 


Integrals  are  constructed  within  a  Galerkin  framework  via 


G?fi 


el 

h  • 


Remark  1:  If  advection  dominates  diffusion,  unless  h  <  v  or  the  ex¬ 
act  solution  is  globally  smooth,  the  Galerkin  approximation  expressed  by 
(2)  is  such  that  qh  will  suffer  from  severe/unacceptable  oscillations  [11,  39]. 
Furthermore,  if  the  discretization  relies  on  high-order  methods,  Gibbs  os¬ 
cillations  may  occur  regardless  of  the  size  of  the  grid.  Different  ways  to 
improve  stability  will  be  described  in  the  following  sections. 


2.1.  The  spectral  element  method 

The  weak  form  in  Eq.  (2)  is  used  with  quadrilateral  elements  of  order 
p,  where  the  element-wise  solution  qh  is  approximated  by  the  expansion 
J2k=i  ^fc(x)  on  Np  =  (p  +  l)2  collocation  points  within  the  element. 
The  expansion  functions  are  constructed  as  the  tensor  product  from  the 
Lagrange  polynomials  /ij(£(x))  and  /ij(r/(x))  of  order  p  as: 

^fe(x)  =  /ti(£(x))  <g>  hj(ri(x.)),  Vi,  j  =  1,  ...,p  +  1.  (4) 

/ij(£(x))  and  hj(ri(x))  are  the  polynomials  associated  with  the  LGL  points 
and  rjj,  respectively.  The  LGL  points  are  the  zeros  of 

(i-eV*(0  =  o 

where  P'N  is  the  derivative  of  the  Nth- order  Legendre  polynomial.  Quadra¬ 
ture  is  performed  on  the  reference  element  Clh  =  [ — 1 , 1] 2  with  LGL  points 
that  have  quadrature  weights  oj.  Substitution  of  the  expansion  J2k= l  ^fc(x)  Qk(t) 
into  the  weak  form  (2)  yields  the  semi-discrete  (in  space)  matrix  problem 

r)nh 

M-^  +  Aqh  +  Dqh  =  0  (5) 

ot 
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where  q/f  is  the  array  of  the  unknowns  on  the  grid  points,  and  M,  A,  and  D 
are  the  global  mass,  advection,  and  diffusion  matrices,  respectively.  These 
matrices  are  obtained  from  the  direct  stiffness  summation  (DSS)  of  the  ele¬ 
mental  matrices  Mei,  Ael,  and  Dfi  given  by: 


m %\=  f  T/,,'1'/  d,nel 

(6a) 

Jci* 

A ek\  =  [  u  •  V'kfc  dnel 

JQel 

(6b) 

D  %[=  f  i A/'Pfc.V'P/dfU 

Jnel 

(6c) 

The  mass  matrix  M  is  diagonal  assuming  inexact  integration. 

All  the  integrals  defined  above  are  approximated  by  the  quadrature  for¬ 
mula 

r  rl  rl  A  ^  ^  A 

n^dx=  /  (7) 

J-1J-1  1=1  j=\ 

where  J  is  the  Jacobian  matrix  associated  with  the  map  between  the  physical 
element  £lh(x,  y )  and  the  reference  element  Clh(£,  v)-  The  integration  is  exact 
up  to  polynomials  of  order  2p-l.  The  p+1  LGL  points  lie  along  the  edges 
and  in  the  interior  of  the  elements.  For  more  on  SEM  see,  e.g.,  [40,  38]. 

We  solve  the  linear  system  of  ordinary  differential  equations  (5)  with  an 
appropriate  strong-stability  preserving  (SSP)  time  integrator.  In  particular, 
we  use  a  five-stage  explicit  third-order  Runge-Rutta  method  (RK35)  [41]. 
SSP  methods  avoid  the  production  of  additional  oscillations  or  damping. 

From  3rd-order  and  up,  the  disposition  of  the  nodes  of  spectral  elements 
differ  from  classical  finite  elements  in  the  way  represented  in  Figure  1  for  a 
4th -order  element. 
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Figure  1:  Nodes  disposition  for  4t,l-order  FE  (left),  and  SE  (right). 


2.2.  Stabilization  techniques 

Aliasing  and  Gibbs  oscillations  render  the  Galerkin  solution  of  (1)  un¬ 
stable  since  unwanted  oscillations  will  pollute  the  numerical  solution.  In 
the  framework  of  spectral  elements,  the  common  strategy  to  control  these 
oscillations  is  the  use  of  anti-aliasing  filters  (See,  e.g.,  [42,  43,  44,  38,  45] 
and  references  therein).  Filtering,  however,  suffers  from  non-positivity  that 
is  unacceptable  in  most  problems  that  involve  transport.  A  suitable  alterna¬ 
tive  to  filters  may  be  the  use  of  a  stabilized  spectral  element  approximation 
based  on  a  residual-based  diffusion-like  term  added  to  the  left  hand-side 
( LHS)  of  (2).  A  stabilization  technique  should  have  the  important  property 
of  consistency  (for  instance,  artificial  diffusion  stabilizes  but  is  not  necessar¬ 
ily  consistent).  Thus,  the  additional  viscous  term  should  vanish  as  the  size 
of  the  element  approaches  zero.  The  stabilized  counterpart  of  (2)  is: 


(^,-jL)+aWh,qh)  +  bW\qh)  =  WhJ)  V 4>h  €  Wh,  (8) 

where  b(i^h,qh)  is  the  stabilization  term. 

A  general  method  for  finding  r  does  not  exist,  and  different  choices  for 
t  dramatically  influence  the  accuracy  of  the  solution.  The  literature  on  the 
optimal  selection  of  r  is  vast,  and  we  refer  to  [25,  26]  for  a  comprehensive 
analysis  of  different  definitions. 
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Method 


AV/HV 

SU 

SUPG 

GLS 

VMS 


Table  1:  Stabilized  methods 

b^'\qh) 

fnh  c^h)v  [vV]  <mh 

Sah  Ctyh)u  [u  ■  Vqh]  <mh 


L  Wh)T 


%-+* -Vqh 


-L  £*Wh)T 


^  +  u-Vqh 


^  +  u  ■  Vg' 


uAqh  -  / 
vAqh  -  / 

-  vA qh  -  f 


dnh 

dflh 

dQh 


£  =  Vaiph 


£  =  u  ■  Vtph 


£  =  u  • 


£  =  u-  X7iph  -  vAiph 


£*  =  — u  •  Viph  -  vA ip 


In  table  1,  v  indicates  a  constant  diffusion  coefficient,  and  a  =  f3/2,  where 
/3  is  a  positive  even  power  of  the  hyper-viscosity  operator  (/3  =  2  yields 
the  usual  AV).  Although  HV  is  scale-selective  (i.e. ,  it  damps  only  higher 
frequencies),  it  is  not  consistent,  nor  is  it  physical.  In  fact,  to  maintain 
the  correct  physical  dimensions  of  the  hyper-viscous  operator,  the  value  of 
the  diffusivity  coefficient  v  must  be  different  when  different  a  are  used.  Its 
selection  is  hence  not  trivial.  Furthermore,  as  Figure  20 -d  indicates,  the 
diffusion  is  isotropic  and  spatially  homogeneous.  The  operator  does  not 
incorporate  the  problem’s  physics. 

The  method  of  Douglas  and  Wang  (DW)  [46]  was  omitted  as  it  can  be 
included  into  the  VMS  method,  although  its  derivation  was  specifically  done 
in  the  context  of  Stokes  problems  rather  than  scalar  advection-diffusion.  The 
SUPG  and  VMS  will  be  described  in  the  following  subsections. 

2.2.1.  Streamline-upwind/ Petrov-  Galerkin  (SUPG) 

The  SUPG  method  was  designed  by  Brooks  and  Hughes  [9]  and  was 
later  generalized  for  multidimensional  problems  by  Hughes  and  Mallet  [47]. 
It  is  a  consistent  alternative  to  the  artificial  diffusion  approach  or  to  the 
overly  diffusive  streamline  upwind  (SU)  method.  Its  use  has  been  ubiquitous 
in  the  solution  of  transport  problems  by  the  finite  element  method  (See, 
e.g.,  [48,  28,  49,  50,  51]).  The  application  of  this  strategy  to  higher-order 
schemes  was  first  tested  for  spectral  methods  by  Canuto  and  coworkers  in 
[17,  18,  19,  52],  and  later  by  Hughes  and  coworkers  in  [21]  using  non-uniform 
rational  B-splines  (NURBS).  In  this  paper  we  show  its  properties  when  used 
with  high-order  spectral  elements.  SUPG  is  a  Petrov-Galerkin  method  in 
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that  it  does  not  assume  that  the  basis  and  test  functions  live  in  the  same 
space.  We  introduce  the  additional  space  of  test  functions  wh  defined  by 

=  \wh  :  wh  =i/jh  +  TU-  Vi/Jh  :  ifh  G  Wh}  . 

We  have  the  problem  of  finding  the  function  qh  G  (Wh;  0,  t)  such  that 


(ifh +tu-V  ifh , 


dqh 

~dt 


)+a(‘iph+TVL-'V'il>h ,  qh) 


(l/jb+TU-Vlpb,  f) 


g  wh. 

(9) 


Rearrangement  of  (9)  yields 


r)ah 

-JL)  +  a(il>h,  qh)  -  tyh,  f)  +b^h,  qh)  =  0  G  Wh  (10) 


Galerkin 


where 


b^h,qh)=  f 

■hi 


^  +  u-Vqh-  uAqh  -  f 


t  u  •  \7iph  dQh  (11) 


is  the  stabilizing  term.  In  (11),  dtqh  +  u  •  V qh  —  uAqh  —  f  is  the  residual  of 
the  governing  equation,  and  r  G  L°°  is  the  stabilization  parameter. 


2.2.2.  The  variational  subgrid  scale  formulation  (VMS) 

Following  [33],  we  define  the  idempotent  linear  projector  V  such  that 
Range(V )  =  W  C  W.  and  such  that  ^  G  W  is  the  projection  of  ijj  G  W 
given  by  if  =  Vif.  If  we  denote  W  as  the  space  of  resolved  scales,  and 
build  the  space  W'  that  completes  W  in  W  and  that  we  will  call  the  space 
of  subgrid  scales,  we  have  W  =  W  ©  W',  where  ©  is  the  overlapping  sum 
decomposition  of  the  two  spaces.  Using  the  properties  of  V,  the  definition 
of  W,  and  the  linear  independence  of  W  and  W' ,  we  decompose  the  weak 
form  (2)  into 

($,  If)  +  (^’  ~§t)  +  °^’  ^  +  a^’  Q' ^  =  ^  (12a) 

(V,  +  a{ip',  q)  +  q')  =  (if',  /).  (12b) 
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We  back-integrate  by  parts  the  bilinear  forms  a(-,  •)  in  (12)  that  depend  on 
the  subgrid  scales,  yielding 

($,  ft )  +  (^’  lit)  +  a^’  ^  ~  In  =  (“05  /)  (13a) 

(V,  ft) + X  ^ + X  ^ = ^  <'i3b') 

where  all  the  boundary  fluxes  and  subscales  vanish  on  the  boundary,  and 
C*  :=  —  —  u- V  —  vAq  is  the  adjoint  operator  of  C.  Equating  the  integrands 

in  (13b),  we  obtain: 


£(q')  =  f-£(q),  (14) 

where  the  time-dependent  contributions  are  included  in  C.  The  solution  of 
(14)  for  q'  must  be  inserted  into  (13a)  to  compute  the  large-scale  solution  q 
by  means  of  the  finite  or  spectral  element  method. 

One  fundamental  step  in  the  derivation  of  the  subgrid  scales  is  the  ap¬ 
proximation  of  q'  from  Equation  (14).  In  the  context  of  high-order  spectral 
elements  with  non-equispaced  nodes,  we  adopt  the  approach  of  Houzeaux, 
Eguzkitza,  and  Vazquez  [34]),  who  used  bubble  functions  to  approximate  q' 
in  the  solution  of  equation  (13b).  To  do  this,  we  use  (14)  and  re-express 
(13b)  using  the  space  of  bubbles  that  vanish  at  the  nodes  of  every  element. 
We  have: 


ifr'  R(q)  d£lh 


W  €  Wq, 


(15) 


where  R(q)  =  dtqh  +  u  •  S7qh  —  uAqh  —  f  is  the  residual  of  the  governing 
equation,  and  Wq  indicates  the  space  of  bubbles  (i.e. ,  subgrid-scale  functions 
that  are  zero  at  the  element  nodes,  and  hence  ensure  nodal  exactness  of  the 
numerical  solution).  We  omitted  the  time-dependent  terms  since  we  used  a 
sufficiently  small  At.  (15)  implies  the  strong  form 


C(q')  =  R(q), 


(16) 


to  be  solved  for  q1  with  Dirichlet  boundary  conditions:  ^(O)  =  0  and 
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q'[h)  =  0,  on  every  element  of  length  h.  Using  the  residual-free  bubble 
b(x)  ([53,  54]),  we  substitute  the  expression  q'(x)  =  b(x)R(qh)  into  (15)  and 
solve  for  b{x)  with  boundary  conditions  6(0)  =  0  and  b{h)  =  0.  We  have: 
C(b(x)R(qh))  dVth  =  R(qh).  To  proceed,  we  assume  that  R(qh)  is  constant 
(i.e.,  we  can  think  that  R(qh)  is  always  known  from  the  previous  time-step), 
and  take  it  out  of  the  operator  £(■).  This  yields  C(b(x))  =  1,  that,  for  the 
one-dimensional  steady-state  advection-diffusion  equation 

ubx(x)  -  vbxx(x)  =  1,  (17) 


has  solution 


_  x  hi-  exu!v 
u~^~  u  ehu/u 


(18) 


For  linear  elements,  h  is  simply  the  length  of  the  element.  For  higher-order 
finite  elements,  h  becomes  a  fraction  of  the  total  element  size  h  if  the  internal 
nodes  are  equispaced.  In  the  case  of  spectral  elements,  where  the  LGL  points 
are  unevenly  distributed,  (18)  is  computed  by  using  h  as  the  local  distance 
between  two  consecutive  points.  The  distributin  of  b(x)  along  the  element  is 
represented  in  Figure  2.  We  use  the  approximation  q'  ~  rR  with  r  defined 
in  (23),  substitute  it  into  (13a)  to  find  the  VMS  (or  generalized)  stabilized 
method 


(^’ ft)  +a^,q^~  In  £*WrR(cl =  (/,# 


(19) 


Eq.  (19)  differs  from  the  weak  form  in  Eq.  (2)  by  the  addition  of  a 
stabilization  term  which  models  subgrid  scales.  The  additional  third  term 
is  the  viscous-like  contribution  that  stabilizes  the  equation. 

Observations  on  time-dependent  subgrid-scales:  The  time-dependent 
approximation  (13)  includes  a  contribution  from  the  time  evolution  of  the 
subscales  given  by  dtq' .  This  may  require  tracking  of  the  subscales  in 
the  time  integration,  unless  the  hypothesis  of  quasi-static  subscales  (i.e. 
dtq'  ~  0)  is  applied  (See  [36]  for  details).  Under  this  hypothesis,  the  con¬ 
tribution  from  the  subgrid  scales  only  appears  in  the  steady  part  of  the 
Galerkin  approximation.  If  a  sufficiently  small  time-step  is  used  with  an 
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explicit  time  integrator,  we  do  not  lose  accuracy  with  the  quasi-static  hy¬ 
pothesis.  With  the  use  of  large  time-steps  with  semi-implicit  time  integrators 
in  atmospheric  simulations,  tracking  of  the  subscales  is  hence  needed.  This 
issue  is  reserved  for  future  work  by  the  authors. 

Remark  3:  It  must  be  pointed  out  that,  for  the  solution  of  the  scalar 
transport  equation,  the  expressions  for  SUPG  and  VMS  are  the  same.  Any 
difference  is  caused  by  a  different  definition  of  r. 

2.3.  The  intrinsic  time  r 

A  standard  intrinsic  time  r  for  finite  elements  was  derived  by  Franca, 
Frey  and  Hughes  [28]: 

r  =  (2°) 

2||U||2 

where  £  is  a  function  of  the  Peclet  number  Pek  which  relates  the  advection 
effects  with  respect  to  diffusion  via: 


d  hk\\u\\2 

Pek  =  mk — — — •  (21) 

In  (21)  mk  is  an  algorithmic  constant  that  should  contain  information  re¬ 
garding  the  interpolating  functions.  For  linear  elements,  one  definition  of 
£  was  derived  using  nodal  exactness  in  [9],  i.e.  the  finite  element  solution 
using  £  defined  in  this  way  is  equivalent  to  the  analytic  solution  at  the  nodes 
of  the  finite  element  mesh  of  a  one  dimension  problem.  In  the  context  of 
subgrid  scale  approximations,  this  condition  is  satisfied  when  the  subgrid 
scales  are  zero  on  the  element  nodes.  We  write: 

£(Pefc)  =  coth(Pefc)  -  — .  (22) 

Pek 

A  whole  set  of  definitions  of  £  appears  in  the  literature:  see  a  review  in  [34]); 
for  systems  of  PDEs,  see  [55]. 

2.3.1.  t  for  spectral  elements 

It  is  assessed  that  for  one-dinrensional  linear  elements,  a  subgrid  scale 
method  gives  nodally  exact  results  if  hk  is  the  full  distance  between  the 
boundary  nodes  of  the  element.  To  derive  it  in  the  context  of  higher-order 
spectral  elements  with  non-equispaced  nodes,  we  extend  the  approach  of 
Houzeaux,  Eguzkitza,  and  Vazquez  [34]  developed  for  quadratic  and  cubic 
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one-dimensional  finite  elements,  using  bubble  functions  in  the  approximation 
of  the  subgrid  scales  to  solve  equation  (13b).  The  stabilization  parameter 
t  is  built  inside  the  element  as  a  function  of  the  bubbles  on  every  segment 
delimited  by  two  consecutive  nodes.  We  write: 


T*+1  = 


T  1)  Jxlgl(i) 


b(x)  dx, 


(23) 


where  xigi(i  +  1)  and  xigi(i )  are  the  coordinates  of  two  consecutive  LGL 
points.  Using  (23),  the  stabilization  integral  is  decomposed  into  a  summa¬ 
tion  of  all  the  contribution  of  the  subelements  [xigi(i  +  1)  —  xigi(i)].  Uniform 
spacing  implies  uniform  t  for  each  sub-element.  The  uneven  spacing  of  the 
element  nodes  is  the  major  difference  with  respect  to  the  definitions  derived 
in  previous  studies.  In  this  case  the  intrinsic  time  is  non-uniform  along  the 
element,  as  it  appears  in  Figure  2,  where  bubbles  and  the  corresponding  r 
are  displayed  for  an  element  of  order  7.  We  expand  r  using  the  Galerkin 
approximation 


p+ i 

r=^2rk'ijjk,  (24) 

k= i 

where  ipk  €  yih  are  the  same  basis  function  previously  introduced  in  the 
approximation  of  the  solution  by  spectral  elements.  The  coefficients  rk  cor¬ 
respond  to  the  value  of  r  at  node  k  and  whose  value  is  obtained  by  (23). 
This  step  allows  the  algorithm  to  maintain  the  simplicity  similar  to  linear 
elements  in  that  we  are  using  one  stabilization  integral  along  the  element; 
however,  the  characteristics  of  the  non-uniform  grid  are  incorporated  into 
the  method.  Throughout  the  paper  we  use  r  in  (23)  for  all  our  simulations. 
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Figure  2:  Bubbles  b(x)  and  r  for  a  7th-order  unitary  spectral  element.  The  biggest  bubble 
in  the  plot  is  the  bubble  that  a  linear  element  would  have. 

2-4-  Spurious  oscillations  at  layers  diminishing  (SOLD)  methods 

Methods  in  the  form  of  (19)  may  produce  overshoots  and  undershoots 
in  the  proximity  of  an  internal  or  boundary  layer.  These  unwanted  oscilla¬ 
tions  can  be  suppressed,  without  affecting  the  global  solution,  by  adding  an 
additional  diffusive  term  of  the  form 

Wfc,f  Vqh),  (25) 


where  consistency  must  be  respected  through  a  proper  construction  of  f. 
We  would  like  to  have  a  method  that  does  not  modify  the  diffusion  in  the 
streamline  direction  since  that  is  already  accounted  for  by  the  stabilization 
term,  but  also  avoids  overdamping  in  the  crosswind  direction.  The  compre¬ 
hensive  set  of  tests  performed  by  John  and  Knobloch  reveals  that  Codina’s 
[27]  is  among  the  best  methods  that  satisfy  these  conditions  when  used  with 
finite  elements.  In  [27],  f  is  defined  by: 


1 

f  =  -max 
2 


_2 v_\  \R(qh)\  (  u0u 
u\\\hk  J  k  I  IVi/'  II  V  |u|2 


(26) 
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where  C is  a  constant,  U||  is  the  velocity  component  in  the  streamline  direc¬ 
tion,  and  <8>  indicates  a  tensor  product.  Codina  suggests  C  =  0.7  for  linear 
and  bilinear  elements,  and  C  =  0.35  for  quadratic  and  biquadratic  elements. 
However,  for  higher  order  elements  using  LGL  points,  we  found  that  the  best 
results  were  obtained  by  setting  C  =  1,  as  long  as  hk  is  selected  properly  in 
the  construction  of  both  f  and  t*,. 

An  alternative  to  (25)  comes  from  Johnson,  Schatz,  and  Wahlbin  [56] 
who  defined  the  following: 

(fux  •  VV^,  ux  •  Vqh)  ,  ux  =  (~W\U\  (27) 

|u| 

In  the  current  work,  (27)  gives  better  results  than  (25),  and  was  then  used 
throughout.  The  results  obtained  with  this  technique  are  labeled  with  DC 
for  Discontinuity  Capturing. 

3.  p- adaptivity 

p-adaptivity  is  one  additional  tool  that  can  further  help  the  suppression 
of  Gibbs  oscillations.  The  concept  is  simple  and  is  easily  coded  on  structured 
grids.  It  consists  of  identifying  the  position  of  discontinuous  solutions,  and 
dropping  the  order  of  discretization  to  1st  for  all  the  elements  that  fall 
within  the  discontinuity.  The  discontinuity  is  sought  with  a  proper  error 
estimator.  The  simple  physics  of  the  advection-diffusion  problems  discussed 
below  allows  for  the  energy-norm  of  the  gradient  of  the  solution  to  be  a 
sufficiently  good  estimator  for  the  current  study.  However,  more  advanced 
methods  should  be  considered.  Algorithm  1  is  a  simple  implementation  of 
this  concept  within  our  code.  The  overhead  is  none  because  the  number  of 
loop  operations  does  not  change  with  respect  to  the  original  case  where  high- 
order  elements  are  used  everywhere.  Clearly,  the  algorithm  can  be  easily 
optimized,  but  we  present  the  pseudo-code  below  for  the  sake  of  clarity. 
The  method  was  applied  to  a  two-dimensional  advection-diffusion  problem 
with  internal  and  boundary  layers  in  a  skew  velocity  field.  Results  are  shown 
in  Figure  9. 

4.  Numerical  testing 

The  algorithms  discussed  in  Section  2  were  tested  using  standard  prob¬ 
lems  described  in  the  literature  of  stabilized  finite  element  methods  for  the 
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Algorithm  1  Compute  the  ls<-order  rhs 

//  Loop  over  all  elements  (iel)  of  the  high-order  computational  grid  : 
for  iel  =  1  to  nelem  do 

/ /  Check  if  the  element  contains  a  discontinuity: 
if  iel  is  s.t.  1 1 V ^ 1 1 2  >  e  then 

/ /  Treat  element  iel  as  a  grid  of  ( ngl  —  1)  x  ( ngl  —  1)  sub-elements: 
for  isubel  =  1  to  maximum  number  of  sub-elements  do 

Create  rhs  using  ls<-0rder  basis  functions  and  integration  rule 

end  for 

else 

//  Do  as  usual: 

Create  rhs  for  the  high-order  spectral  element. 

end  if 

end  for 


advection-diffusion  (AD)  equation: 

•  ID  Steady-state  homogeneous  advection-diffusion  (St- ID) 

•  ID  Steady-state  advection-diffusion  with  source  ( St-ID-S) 

•  2D  Steady-state  advection-diffusion  with  internal  and  boundary  layers 

(St- 2D). 

•  2D  Time-dependent  advection-diffusion  with  “L”-shaped  discontinuity 

(Trl-2D). 

•  2D  Time-dependent  advection  of  a  sharp  tracer  in  a  doubly  periodic 
channel  (Tr2-2D). 

St-1D:.  One-dimensional  steady-state  advection-diffusion.  The  tracer  qh  is 
propagated  with  constant  velocity  u  =  1m  s_1  and  diffusivity  u  =  1/512  m 2  s-1 
first  on  two  elements  of  order  p  =  10  (Figure  3),  and  then  on  four  elements 
of  order  p  =  12  (Figure  4).  The  domain  is  the  line  segment  D  =  [—1,1]  with 
Dirichlet  boundary  conditions  qh(— 1)  =  0  and  g^(l)  =  1.  We  compared 
the  filtered  (top  row)  against  the  stabilized  solution  (bottom  row)  and  ob¬ 
serve  a  decrease  of  oscillations  and  undershoots.  Also,  at  higher  order  and 
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finer  resolution,  the  capabilities  of  the  filter  are  clearly  being  challenged 
by  the  presence  of  the  boundary  layer  at  x  =  1.  At  the  same  time,  small 
oscillations  near  the  nodes  of  the  element  by  the  boundary  layer  are  not 
completely  suppressed  by  the  stabilized  method  either;  hence,  additional 
localized  smoothing  is  sought. 


(a)  Filter.  2  el,  p  =  10 


(b)  VMS.  2  el,  p  =  10 

Figure  3:  St-ID:  v  =  l/512m2s-1.  The  ex¬ 
act  solution  is  dashed.  The  circles  indicate 
the  grid  points. 


(a)  Filter.  4  el,  p  =  12 


(b)  VMS.  4  el,  p  =  12 

Figure  4:  St-ID:  v  =  1/512  m2  s-1.  The  ex¬ 
act  solution  is  dashed.  The  circles  indicate 
the  grid  points. 


St-1D-S:.  Steady  state  advection-diffusion  with  source  term  /  =  1.  qh  is 
propagated  with  constant  velocity  u  =  lm  s”1  and  two  different  diffusivities: 
v  =  5  x  10”3  m 2  and  v  =  5  x  10~2  rn2  s  .  The  domain  is  the  line  segment 
17  =  [0, 1]  and  homogeneous  Dirichlet  boundary  conditions  are  imposed. 
The  domain  is  subdivided  into  two  elements  of  order  p  =  16  and  runs  are 
compared  using  filtered  SE  (top  row  in  Figures  5  and  6),  and  VMS  (bottom 
row).  We  observe  a  very  similar  behavior  of  the  solution  among  the  two 
different  cases  in  the  smooth  problem  (Figure  5).  The  results  are  comparable 
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to  the  ones  obtained  by  Houzeaux  et  al.  with  their  r  for  quadratic  and  cubic 
elements  in  [34]. 


(b)  VMS  (b)  VMS 


Figure  5:  St-ID-S:  v  =  5  x  1CP3  m2s_1.  2 
16t/l-order  elements.  The  exact  solution  is 
dashed.  The  circles  indicate  the  grid  points. 


Figure  6:  St-1D-S\  v  =  5  X  10-2  m2s_1.  2 
16t,l-order  elements.  The  exact  solution  is 
dashed.  The  circles  indicate  the  grid  points. 


St- 2D:.  Standard  steady  advection-diffusion  skew  to  the  mesh  (e.g.,  [27]): 
a  discontinuity  is  propagated  with  constant  velocity  u  =  (l,-2)mV1  and 
diffusivity  v  =  1CP8  m2  s-1  in  the  unit  square  Q  =  [0, 1]  x  [0, 1].  The  initial 
configuration  is  shown  in  Figure  7. 
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Figure  7:  St-2D:  initial  configuration  of  the  steady-state  problem 
Dirichlet  boundary  conditions  are  prescribed: 

fl  if  y  =  1, 

Qh  =  |  1  if  x  =  0  and  y  >  0.7, 

I  0  otherwise. 

In  Figures  8-12  we  illustrate  the  run  of  the  same  case  with  different  number 
of  elements  and  order  of  the  interpolating  polynomials.  For  direct  compari¬ 
son  of  our  solution  with  the  ones  in  the  existing  literature  of  finite  elements, 
we  first  run  the  test  with  linear  elements  (p  =  1),  and  present  the  results  in 
Figure  8.  The  multiscale  solution  of  this  problem  (see  Figure  8a)  shows  im¬ 
portant  boundary  and  internal  layers  that  are  damped  with  the  discontinu¬ 
ity  capturing  techniques  of  Section  2.4.  The  application  of  the  discontinuity 
capturing  (DC)  scheme  greatly  improves  the  solution  and  yields  monotonic¬ 
ity  (see  Figure  86).  In  Figure  9  we  maintained  the  same  number  of  nodes 
of  the  previous  run,  but  increased  to  4th  the  order  of  interpolation  to  as¬ 
sess  the  algorithm  in  the  context  of  this  paper  (i.e.  50  elements  of  order 
4  were  used  instead  of  200  elements  of  order  1).  The  similar  behavior  of 
the  solution  with  respect  to  the  lst-order  polynomial  run  suggests  that  the 
residual-based  methods  as  implemented  in  this  study  may  not  be  sensitive  to 
the  distribution  of  the  interpolation  nodes  within  the  elements  edges.  As  it 
appears  in  Figures  9  c,  the  behavior  is  completely  analogous  to  the  previous 
run.  However,  monotonicity  is  lost  in  two  singular  nodes:  with  reference 
to  Figure  9c,  the  4^'-order  solution  is  smooth  and  monotone  everywhere 
except  for  the  nodes  represented  as  points  A  and  B  in  Figure  7.  This  is 
not  surprising:  at  A  and  B  the  tracer  is  leaving  the  boundary  with  a  skew 
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angle;  an  incorrect  imposition  of  boundary  conditions  at  these  nodes  may 
be  causing  the  problem.  The  numerical  singularity  at  this  points  should 
be  addressed  but  it  will  not  be  done  in  the  current  work.  These  are  fully 
suppressed  by  applying  the  p-adaptivity  algorithm  described  above,  as  it  is 
shown  in  Figure  9  d. 

Decreasing  the  number  of  computational  nodes  by  doubling  the  order 
from  4  to  8  and  setting  the  number  of  elements  to  10  in  x  and  z,  even  with  a 
discontinuity  capturing  term,  the  solution  starts  to  lose  monotonicity.  This 
appears  in  Figures  11  and  12,  where  extrema  get  larger  than  in  the  previous 
cases.  This  problem  shows  that  the  construction  of  the  stabilizing  parameter 
r  should  include  information  on  the  order  of  the  interpolating  polynomial. 

For  a  better  view  of  the  problem,  in  Figures  10  and  12  we  present  a  ver¬ 
tical  slice  of  the  solution.  The  boundary  layers  are  evident.  Their  damping, 
however,  is  clear  if  V MS  +  DC,  and  possibly  p-adaptivity,  are  applied. 


(a)  VMS  (b)  VMS  +  DC 

Figure  8:  St-2D :  steady-state  solution  on  200x200  lst-order  elements.  (For  plotting  only, 
the  data  are  interpolated  to  a  50  x  50  node  grid  using  Octave  [57]). 
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(Mm,  Max)  =  (-5,509e-3, 1.026)  (Min,  Max)  =  (-2.161e-07, 1) 

(c)  VMS  +  DC  (d)  VMS  +  DC  +  p-adaptivity 

Figure  9:  St-2D :  steady-state  solution  on  50x50  4th-order  elements. 


Figure  10:  St-2D\  steady-state  solution  on  50x50  4t,l-order  elements.  Vertical  slice  at 
2  =  0.3 
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(Min,  Max)  =  (-3.307, 4.809) 


(Min,  Max)  =  (-2.456e- 1 ,  1.368) 


(a)  Filter 


(b)  VMS 


(Min,  Max)  =  (-7.386e-3,  1.016) 


(Min,  Max)  =  (-3.724e-7,  1.005) 


(c)  VMS  +  DC 


(d)  VMS  +  DC  +  p-adaptivity 


Figure  11:  St-2D :  steady-state  solution  on  10x10  8{,1-order  elements. 
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Figure  12:  St-2D\  steady-state  solution  on  10x10  8t,l-order  elements.  Vertical  slice  at 
2  =  0.3 

Trl-2D:.  Transient  advection-diffusion  of  an  L-shaped  discontinuity  in  a 
flow  where  v  =  10~6  m 2  s_1  and  the  velocity  u  of  magnitude  |u|  =  0.5\/2  m  s~ 
is  at  45°  with  respect  to  the  axis  ( x,z ).  The  initial  configuration  is  shown 
in  Figure  13. 


qh  =  i  b  qh  =  o 


Figure  13:  Trl-2D\  initial  configuration  of  the  L-shaped  problem 

The  convex  shape  of  the  sharp  discontinuity  makes  this  problem  more 
challenging  than  the  previous  case  [58],  and  is  chosen  to  analyze  robustness 
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and  accuracy  of  the  algorithm.  Runs  were  performed  at  two  different  reso¬ 
lutions  and  two  different  orders  of  interpolating  polynomials.  In  particular 
we  have:  approx.  100  points  per  side  using  25x25  4 th  —  order  elements  (Fig¬ 
ure  14),  and  12x12  84,l-order  elements  (Figure  15);  and  approx.  200  points 
per  side  using  50x50  4<fc-order  (Figure  17),  and  25x25  8<ft-order  (Figure 
17).  In  the  figures,  Filter  means  that  the  SEM  solution  was  filtered  at  every 
time-step.  VMS  and/or  DC  indicate  that  the  SEM  solution  is  stabilized  by 
the  VMS  with  or  without  a  discontinuity  capturing  term  (DC).  VMS  +  DC 
+  p-adaptivity  indicates  the  contribution  of  p-adaptivity  as  well.  Positiv¬ 
ity  is  not  preserved  in  the  solution  obtained  with  a  filter.  The  sharp  front, 
in  fact,  makes  the  filter  inappropriate.  However,  similarly  to  the  steady 
advection-diffusion  test  St-2D,  the  VMS-stabilized  solution  of  this  problem 
is  characterized  as  well  by  the  formation  of  internal  layers  that  run  along  the 
edges  of  the  tracer  in  the  direction  of  the  flow  (See,  e.g.,  Figure  14b),  and 
VMS  is  not  sufficient  to  preserve  nronotonicity  unless  it  is  supplemented  by 
the  additional  DC  term  defined  in  (27).  This  effect  is  displayed  in  Figures 
14,15,16,  and  17. 

The  consideration  made  for  problem  St-2D  on  the  singular  peaks  that 
form  at  the  nodes  where  the  tracer  leaves  the  boundary  at  an  angle,  applies 
here  at  nodes  A  and  B  of  Figure  13.  This  is  visible  in  Figure  18  obtained 
by  slicing  the  tracer  along  z  =  0  in  Figures  16  and  17,  respectively.  The 
problem  is  solved  by  the  application  of  p-adaptivity. 

As  the  order  of  interpolation  is  increased  from  4th  to  8th,  the  smooth 
solution  begins  to  lose  positivity.  As  interpreted  for  St-2D,  the  solution  is 
clearly  being  affected  when  the  interpolation  nodes  are  densely  clustered 
towards  the  boundaries  of  the  elements,  as  is  the  case  for  higher  order. 
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rngmum 


(a)  Filter 


(a)  Filter 


(d)  VMS  +  DC  +  p-adaptivity 


(d)  VMS  +  DC  +  p-adaptivity 


Figure  16:  Trl-2D :  4' 


order  50x50.  t 


order  25x25.  t 


(a)  4th  —  order  50x50 


(b)  8th  -  order  25x25 


Figure  18:  Trl-2D\  Vertical  slice  at  z  =  0.0. 

Tr2-2D:.  Linear  advection  of  a  2D  square  wave  along  x  in  the  periodic 
domain  D  =  [0,1]  x  [0,1]:  the  tracer  is  transported  with  velocity  u  = 
(1/2,  0)  ms-1  for  one  periodic  revolution  along  x.  The  initial  concentration 
qh  =  1  is  centered  at  (xc,zc)  =  (0.5,  0.5)  (Figure  19).  The  computational 
finite  domain  consists  of  11  x  11  quadrilaterals  of  order  11. 


Figure  19:  Tr2-2D\  initial  configuration  of  the  pure  advection  problem. 

As  in  the  steady  case,  Figures  20-22  display  improvement  of  the  solution 
in  terms  of  nronotonicity  when  the  VMS  method  is  used  instead  of  the 
filter.  The  combination  of  VMS  and  filtering  is  not  recommended  (result 
not  shown);  although  VMS  alone  controls  the  over-  and  under-shootings 
along  the  streamlines,  the  addition  of  the  filter  at  the  end  of  every  time  step 
degrades  positivity  in  the  neighborhood  of  large  gradients. 

In  Figures  21  and  22  we  present  the  streamline  and  crosswind  sections  of 


the  solution  obtained  by  slicing  the  tracer  along  z  =  0.5  and  x  =  0.5,  respec¬ 
tively.  Unlike  the  previous  problems  characterized  by  internal  and  boundary 
layers,  for  pure  advection  the  VMS  preserves  the  maximum  and  minimum 
concentrations  q^ax  =  1>  and  Qmin  =  0  and  is  free  of  spurious  oscillations. 
As  a  term  of  comparison,  we  present  the  result  of  classical  artificial- viscosity 
in  Figures  20-22 d.  The  exact  solution  of  this  test  at  the  final  time  t  =  2 
corresponds  to  the  initial  condition.  The  normalized  L\,  L2  and  L0 0  errors 
of  the  computed  solution  with  respect  to  the  exact  solution  were  computed 
on  different  grids  to  obtain  the  convergence  curves  for  different  orders  of  in¬ 
terpolation.  Figure  23  shows  a  log-log  plot  of  the  L error  against  the  grid 
size  A(x,  z)  =  0.2, 0.1,  0.05,  0.025, 0.0125,  0.00625  ( nel  =  5x5, ...,  160  x  160), 
where  each  curve  represents  the  error  obtained  with  different  orders  of  in¬ 
terpolation  (p  =  4, ...,  10).  At  a  given  p,  the  At  was  decreased  to  maintain  a 
constant  Courant  number  for  all  resolutions.  At  a  given  p.  the  slope  of  the 
curves  increases  at  increasing  resolution.  This  is  expected  by  the  spectral 
element  solution  in  that  the  error  reduces  exponentially.  However,  as  p  is 
increased,  the  slope  tends  to  straighten  (See  Figure  24),  indicating  that  the 
larger  number  of  nodes  due  to  the  increased  p  contributes  to  setting  the 
convergence  rate  of  the  method  in  the  same  extent  as  high  /^-resolution  does 
when  p  is  lower.  The  data  used  to  generate  the  plots  are  reported  in  Table 
2. 
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(Min,  Max)  =  (-0.1 191 ,  1.11) 


(Min,  Max)  =  (-0.1845,  1) 


(a)  Galerkin 


(b)  Filter 


(Min,  Max)  =  (-4.125e-06, 1) 

(c)  VMS 

Figure  20:  Tr2-2D:  Surface  plot  of  the  concentration  field:  At  =  0.001s  (except  for  HV: 
At  =  0.0002s),  11x11  elements  with  11th  order  polynomials.  Results  at  t  =  2.0  s  (after 
1  periodic  revolution  along  x). 


(Min.Max)  =  (0.0,  0.9482) 
(d)  AD/HV  v  =  0.001  m2  s”1 
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(a)  Galerkin  (b)  Filter 


(c)  VMS  (d)  AD/HV  v  =  0.001  m2  s”1 

Figure  21:  Tr2-2D:  Streamline  cut  at  0.5  m  in  the  y-direction.  At  =  0.001s,  11x11 
elements  with  11th  order  polynomials.  Results  at  t  =  2  s  (after  1  periodic  revolution 
along  x).  Solid  line  indicates  the  computed  solution.  The  dashed  line  is  the  analytic 
solution. 
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(a)  Galerkin  (b)  Filter 


1  - 

0.8 

0.6 

0.4 

0.2 

0  -  - 

0  0.2  0.4  0.6  0.8  1 

y  -  crosswind  direction 


(c)  VMS 


(d)  AD/HV  v  =  0.001  m2  s”1 


Figure  22:  Tr2-2D:  Crosswind  cut  at  0.5  m  in  the  x-direction.  At  =  0.001s,  11x11 
elements  with  11th  order  polynomials.  Results  at  t  =  2  s  (after  1  periodic  revolution 
along  x).  Solid  line  indicates  the  computed  solution.  The  dashed  line  is  the  analytic 
solution. 
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h 


Figure  23:  Tr2-2D:  Log-log  plot  of  the  L0 0  norm  error  vs.  h  for  polynomials  p  4,6,8,10. 
Original  tabled  data  are  reported  in  Appendix  B. 
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Table  2:  Tr2-2D\  Maximum  and  minimum  values  of  concentration  and  errors  in  different 
norms  after  1  full  periodic  revolution  ( t  =  2s).  At  =  0.001s,  11x11  elements  of  11th 
order. 


Run 

Qmin 

Qmax 

Li 

l2 

Loo 

lioomin 

FILTER 

VMS 

-0.1951E+00 

-0.7010E-06 

1.1010 

1.0000 

0.13969 

0.16581 

0.1732 

0.2379 

0.10134E+00 

0.10965E-05 

0.10134E+00 

0.10965E-05 

h 


Convergence  for  p=6 


Convergence  for  p=10 


Figure  24:  Tr2-2D\  Log-log  plots  of  the  Li,L2,  and  Lx  norm  errors  vs.  h.  The  original 
data  are  reported  in  Appendix  B. 


Remark  on  the  use  of  filters  in  the  previous  results:  In  the  current 
work,  filtering  was  applied  in  the  usual  way  that  has  been  used  previously  in 
SE  models  (see,  e.g.,  [44,  59,  60,  61]).  That  is,  the  filtering  coefficients  were 
defined  at  the  beginning  of  the  simulation  and  applied  after  every  time-step 
using  the  same  filter  matrix  for  all  elements.  It  may  be  possible  to  obtain 
better  results  with  filters  if  they  are  constructed  in  a  specific  way  (e.g.,  each 
element  uses  a  different  filter  matrix  that  is  constructed  dynamically)  but 
a  clear  approach  on  how  to  do  this  remains  an  open  topic  since  this  can  be 
viewed  as  a  classical  limiter  but  for  Spectral  Elements  (see,  e.g.,  [62]). 
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5.  Discussion  and  conclusions 


5.1.  Conclusions 

In  this  paper,  we  proposed  the  use  of  the  Variational  Multiscale  method 
(VMS)  to  stabilize  advection-dominated  problems  solved  with  spectral  el¬ 
ements.  The  stabilization  parameter  r  that  appears  in  the  VMS  scheme 
was  computed  to  include  the  characteristics  of  high-order  spectral  elements 
with  (non-equispaced)  LGL  nodes.  Numerically,  we  demonstrated  that  this 
approach  is  a  possible  alternative  to  the  standard  filters  used  in  the  sta¬ 
bilization  of  the  spectral  element  solvers  if  the  suppression  of  unwanted 
under-  and  over-shoots  is  the  main  concern.  In  the  presence  of  internal  and 
boundary  layers,  VMS  was  coupled  with  appropriate  discontinuity-capturing 
techniques  to  damp  any  oscillations  in  the  proximity  of  such  discontinuities. 
Stabilization  by  these  methods  is  obtained  by  introducing  a  diffusion-like 
term  that,  unlike  hyper- viscosity,  only  acts  in  the  regions  where  oscilla¬ 
tions  occur  (i.e.  large  gradients);  the  solution  is  not  affected  in  smooth 
regions.  Where  needed,  the  combined  action  of  VMS  and  polynomial  adap¬ 
tivity  yields  encouraging  results  for  high-order  spectral  elements.  The  al¬ 
gorithms  were  evaluated  on  a  set  of  standard  tests  of  increasing  difficulty. 
A  significant  improvement  was  observed  in  the  performance  of  the  spectral 
element  solver  as  far  as  the  control  of  maxima  and  minima  is  concerned, 
both  in  the  purely  advective  and  in  the  advective-diffusive  regimes.  The 
most  important  features  of  this  new  approach  are: 

•  Unlike  hyper-viscosity,  the  subgrid-scale  diffusion  only  affects  regions 
of  the  ffow  where  stabilization  is  required. 

•  Under-  and  over-shoots  are  greatly  suppressed  relative  to  traditional 
filters. 

•  The  method  does  not  depend  on  a  free-parameter  assigned  by  the  user. 

5.2.  Application  to  atmospheric  modeling  in  climate  and  weather  prediction, 
and  future  work 

In  [63],  we  implemented  a  Kessler  microphysics  scheme  [1]  within  a  spec¬ 
tral  element  framework,  that  requires  the  advection  of  three  moisture  vari¬ 
ables  (vapor,  cloud,  and  rain  mixing  ratios).  This  microphysics  scheme  will 
be  implemented  in  our  Nonhydrostatic  Unified  Model  for  the  Atmosphere 
(NUMA)  [64]  in  order  to  simulate  both  mesoscale  and  synoptic-scale  at¬ 
mospheric  phenomena.  As  is  well-known,  Galerkin-based  methods  yield 
1)  higher-order  accuracy  and  2)  excellent  dispersion  properties,  which  are 
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both  desirable  for  advection  schemes;  however,  the  resulting  Gibbs  oscil¬ 
lations  produce  strong  gradients  that  must  be  remedied  in  some  fashion. 
At  present,  a  simple-minded  “fixer”  is  applied  whereby  negative  values  of 
the  moisture  variables  are  set  equal  to  zero.  This  fixer  acts  as  an  effective 
mass  source,  thus  violating  the  conservation  properties  of  the  model.  In 
addition,  this  fixer  violates  the  function  space  that  the  spectral  element  so¬ 
lution  inhabits.  For  these  reasons,  monotonic  advection  of  tracer  variables 
is  essential  for  any  atmospheric  model.  The  proposed  VMS  scheme  is  an  ex¬ 
cellent  candidate  since  it  1)  preserves  nronotonicity  better  than  the  standard 
filter  approach,  2)  does  not  significantly  increase  the  cost  of  the  spatial  dis¬ 
cretization  scheme,  and  3)  is  completely  local  in  nature  (i.e. ,  no  additional 
communications  are  required  in  a  parallel  environment),  which  is  necessary 
for  scaling  on  modern  distributed  and  hierarchical  memory  environments. 
In  future  work,  we  will  report  on  the  application  of  VMS  to  NUMA  using 
realistic  mesoscale  test  cases. 
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Appendix  A.  Expression  for  r 

In  this  appendix  we  explicitly  derive  the  expression  for  r  defined  between 
two  consecutive  LGL  points  [xigi(i),  xigi(i  +  1)].  The  bubble  obtained  from 
the  integration  of  (17)  with  boundary  conditions  b(xigi(i))  =  0  and  b(xigi(i  + 
1))  =  0  has  expression: 

x  x(i  +  1)  —  x(i)  ux/v  x{i)eux^l+1^v  —  x(i  +  l)eua:W/i' 

^  ^gux(i+l)  /  v  g  ux{i)/v^  ^guRi+l)  /  v  g  ux(i)/u^ 

The  evaluation  of  the  integral  (23) 
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rf+1  = 


r^Z9z(*+l) 


xigi(i  +  1)  -xigi(i)  Jxlgl{i) 


b(x)  dx 


yields  the  expression: 


x(i+l)  _ 
Tx{i) 


x(i  +  1)  —  x(i) 


x(i)  —  x{i  +  1)  /is  x(i  +  l)  —  x(i)\  eux^+1^l'(x(i)  —  x(i  +  l))2 

u  \u  +  2 


g ux(i)/iy  _  qUx(i-\-1)/u 


When  x(i)  =  0  and  x(i  +  1)  =  h,  we  have  that 

^h_  v  h  heuh//u 
T°  ~  u2  2  u  +  euh/v  -  1  ’ 

from  which,  after  little  algebra,  expression  (20)  is  recovered: 


T  =  A  (coth(pe„)  -  A_)  . 


Appendix  B.  Error  tables 

In  this  Appendix  we  report  the  data  of  error  estimates  from  which  the 
curves  of  Figures  23  and  24  were  generated. 


Table  B.3:  Tr2-2D:  Convergence  results  for  the  transport  problem  using  the  VMS:  p=^. 
At  =  0.001  for  all  the  runs. 


Nelements 

Li 

l2 

Loo 

Loo  min 

Qmax 

Qmin 

5x5 

0.30786 

0.2583 

0.9662E-01 

0.1209E+00 

1.0970 

-0.6290E-01 

10  x  10 

0.19662 

0.2127 

0.4701E-01 

0.4770E-01 

1.0430 

-0.2850E-01 

20  x  20 

0.13288 

0.1845 

0.8739E-02 

0.8737E-02 

1.0090 

-0.6960E-02 

40  x  40 

0.87334E-01 

0.1546 

0.1714E-02 

0.1714E-02 

1.0020 

-0.1793E-02 

80  x  80 

0.61619E-01 

0.1323 

0.4408E-04 

0.4408E-04 

1.0000 

-0.4756E-04 

160  x  160 

0.44715E-01 

0.1146 

0.2650E-07 

0.2649E-07 

1.0000 

-0.4116E-07 
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Table  B.4:  Tr2-2D :  Convergence  results  for  the  transport  problem  using  the  VMS:  p=4- 
Variable  At  for  all  the  runs:  CFL~0. 01025 _ 


Nelements 

At 

Ti 

l2 

Loo 

Loo min 

Qmax 

Qmin 

Mloss 

5x5 

0.0100E-01 

0.3079 

0.2583 

0.9662E-01 

0.1209E+00 

1.0970 

-0.6290E-01 

0.180E-12 

10  x  10 

0.0500E-02 

0.1966 

0.2127 

0.4701E-01 

0.4770E-01 

1.0430 

-0.2850E-01 

0.182E-12 

20  x  20 

0.0250E-02 

0.1329 

0.1845 

0.8739E-02 

0.8737E-02 

1.0090 

-0.6960E-02 

0.350E-12 

40  x  40 

0.1250E-03 

0.8733E-01 

0.1546 

0.1714E-02 

0.1714E-02 

1.0020 

-0.1793E-02 

0.730E-12 

80  x  80 

0.6250E-03 

0.6162E-01 

0.1323 

0.4408E-04 

0.4408E-04 

1.0000 

-0.4756E-04 

0.147E-11 

160  x  160 

0.3125E-04 

0.4470E-01 

0.1136 

0.3661E-07 

0.3661E-07 

1.0000 

-0.5370E-07 

0.603E-11 

Table  B.5:  Tr2-2D:  Convergence  results  for  the  transport  problem  using  the  VMS:  p= 
Variable  At  for  all  the  runs:  CFL~0.0208 

-6. 

Nelements 

At 

£i 

l2 

Loo 

Loomin 

Qmax 

Qmin 

Mioss 

5x5 

0.0100E-01 

0.2186 

0.2348 

0.9349E-02 

0.9349E-02 

1.0090 

-0.8514E-02 

0.9084E-12 

10  x  10 

0.0500E-02 

0.1159 

0.1662 

0.3126E-02 

0.3126E-02 

1.0030 

-0.2680E-02 

0.1810E-11 

20  x  20 

0.0250E-02 

0.1098 

0.1810 

0.8451E-03 

0.8451E-03 

1.0010 

-0.5194E-03 

0.3642E-11 

40  x  40 

0.1250E-03 

0.6983E-01 

0.1357 

0.7971E-04 

0.7971E-04 

1.0000 

-0.9405E-04 

0.7241E-11 

80  x  80 

0.6250E-03 

0.5206E-01 

0.1246 

0.2599E-05 

0.2599E-05 

1.0000 

-0.2007E-05 

0.1454E-10 

160  x  160 

0.3125E-04 

0.3630E-01 

0.1018 

0.1276E-08 

0.1276E-08 

1.0000 

-0.1172E-08 

0.2904E-10 

Table  B.6:  Tr2-2D :  Convergence  results  for  the  transport  problem  using  the  VMS:  p=8. 
Variable  At  for  all  the  runs:  CFL~0.0352 _ 


Nelements 

At 

In 

l2 

Loo 

Qmax 

Qmin 

Mioss 

5x5 

0.0100E-01 

0.1415 

0.18180 

0.2970E-02 

0.2970E-02 

1.0030 

-0.7465E-02 

0.9053E-12 

10  x  10 

0.0500E-02 

0.9237E-01 

0.15087 

0.1852E-02 

0.1852E-02 

1.0020 

-0.1610E-02 

0.1816E-11 

20  x  20 

0.0250E-02 

0.9162E-01 

0.15879 

0.7905E-04 

0.7905E-04 

1.0000 

-0.1115E-03 

0.3612E-11 

40  x  40 

0.1250E-03 

0.5929E-01 

0.12320 

0.2979E-05 

0.2979E-05 

1.0000 

-0.4298E-05 

0.7325E-11 

80  x  80 

0.6250E-03 

0.4432E-01 

0.11340 

0.4332E-07 

0.4332E-07 

1.0000 

-0.1031E-06 

0.1436E-10 

160  x  160 

0.3125E-04 

0.3129E-01 

0.9399E-01 

0.1418E-09 

0.1418E-09 

1.0000 

-0.7009E-10 

0.2854E-10 

Table  B.7:  Tr2-2D :  Convergence  results  for  the  transport  problem  using  the  VMS: 
Variable  At  for  all  the  runs:  CFL~0.0536 

p=10. 

Nelements 

At 

Lx 

l2 

Loo 

Loomin 

Qmax 

Qmin 

Mioss 

5x5 

0.0100E00 

0.1647 

0.2034 

0.3855E-02 

0.3855E-02 

1.0040 

-0.5317E-02 

0.9069E-12 

10  x  10 

0.0500E-01 

0.7786E-01 

0.1415 

0.2796E-03 

0.2796E-04 

1.0000 

-0.2022E-03 

0.1820E-11 

20  x  20 

0.0250E-01 

0.8875E-01 

0.1622 

0.7076E-05 

0.7076E-05 

1.0000 

-0.1081E-04 

0.3651E-11 

40  x  40 

0.1250E-02 

0.5183E-01 

0.1142 

0.1937E-06 

0.1937E-06 

1.0000 

-0.7324E-07 

0.7298E-11 

80  x  80 

0.6250E-03 

0.4083E-01 

0.1115 

0.3612E-08 

0.3612E-08 

1.0000 

-0.1089E-07 

0.1455E-10 

160  x  160 

0.3125E-04 

0.2785E-01 

0.8815E-01 

0.4936E-10 

0.4936E-10 

1.0000 

-0.1127E-10 

0.2842E-10 
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